Robustness of a perturbed topological phase 
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We investigate the stability of the topological phase of the toric code model in the presence of a 
uniform magnetic field by means of variational and high-order series expansion approaches. We find 
that when this perturbation is strong enough, the system undergoes a topological phase transition 
whose first- or second-order nature depends on the field orientation. When this transition is of 
second order, it is in the Ising universality class except for a special line on which the critical 
exponent driving the closure of the gap varies continuously, unveiling a new topological universality 
class. 
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Introduction — The concept of topological quantum 
order was introduced by Wen in the late 1980s, to char- 
acterize the chiral spin state supposed to be relevant for 
high-temperature superconductivity [ ]. Since then, it 
has been shown to be crucial for characterizing differ- 
ent states of matter, among which are fractional quan- 
tum Hall states, and it has become the cornerstone of 
topological quantum computation [2, 3]. Topologically 
ordered quantum systems are mainly characterized by 
a ground-state degeneracy which depends on the Eulcr- 
Poincare characteristic. For connected orientable sur- 
faces, this number is directly related to the genus. Topo- 
logically ordered states cannot be characterized by local 
order parameters and thus fail to be described by Landau 
symmetry-breaking theory. Importantly, this nonlocality 
often implies anyonic statistics and a robustness of the 
corresponding system with respect to any local pertur- 
bation [2, 4, 5], so that they might be used as reliable 
quantum memories [ ] . However, it has early been real- 
ized in the seminal paper of Kitaev [ ] that "Of course, 
the perturbation should be small enough, or else a phase 
transition may occur. " 

The main motivation of the present work is precisely to 
investigate this robustness in the simplest model display- 
ing topological quantum order, namely, the toric code [2], 
and in the presence of the simplest local perturbation, 
i.e., a uniform magnetic field. This model, which might 
be implemented in Josephson junction arrays [7], may 
indeed be considered as the "Ising model of topological 
quantum phase transitions" and has already been studied 
for special directions of the field [8-12] (see also Ref. 13 
for a related problem in Wen's model [14]). Here, we 
address this problem for an arbitrary field direction and 
determine the extension of the topological phase origi- 
nating from the zero-field limit. To compute this phase 
diagram, one faces several difficulties since (i) the lack 
of a local order parameter prohibits any field-theoretical 



approach to analyze the critical properties and (u) one 
can neither perform Monte-Carlo simulations (sign prob- 
lem) nor reliable exact diagonalizations (only small sizes 
are available). Consequently, we combine two different 
techniques. First, we perform high-order series expan- 
sion in the small-field limit using perturbative continu- 
ous unitary transformations (PCUT) [15] and compute 
the ground-state energy as well as the low-energy gap. 
Unfortunately, although such an expansion is very effi- 
cient to characterize second-order transitions [11], it can- 
not locate first-order transitions except in very special 
situations [12]. Second, we use a variational approach 
based on infinite projected entangled pair states (iPEPS) 
[16-18] which is, by contrast, especially sensitive to first- 
order transitions (see, for instance, Ref. 19). Combining 
these two methods, we determined the boundaries of the 
topological phase of the toric code model in an arbitrary 
uniform magnetic field. The resulting phase diagram dis- 
plays many interesting features since, depending on the 
direction of the field, the breakdown of the topological 
phase may be achieved through a first- or a second-order 
transition. In the latter case, the universality class is 
always of Ising type except on a special line where the 
critical exponent driving the closure of the gap varies 
continuously. 

Model and limiting cases — The Hamiltonian of the 
toric code in a uniform magnetic field reads 

s p i 

where A s = Uies^f and B v = Uie P a i (°f' s are the 
usual Pauli matrices). Subscript s (p) refer to sites (pla- 
quettes) of a square lattice and i runs over all bonds 
where spins are located [ ]. Without loss of generality, 
we restrict our study to h a ^ 0, the spectrum being un- 
changed under the transformation h a — > —h a . 

In the zero- field limit, H is exactly solvable since 
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[A s , Bp] = 0. As shown in Ref. 2, the ground-state degen- 
eracy depends on the surface topology so that the system 
is topologically ordered. In this limit, the ground-state 
energy per spin is eo = — J. Elementary excitations are 
obtained by acting onto the ground states with of (charge 
excitations) or erf (flux excitations) operators which lo- 
cally change the eigenvalues of A s or B p . On a torus, only 
pairs of such elementary excitations can be created so 
that, in this case, one has an equidistant spectrum with 
an energy gap A = 4J. By contrast, for open bound- 
ary conditions, the gap is A = 2J since one can create 
states with only one charge or only one flux. Charges 
and fluxes behave individually as hard-core bosons but 
have mutual anyonic (semionic) statistics [2]. In the op- 
posite limit J = 0, the ground state is unique and fully 
polarized in the field direction whatever the boundary 
conditions so it is not topologically ordered anymore. It 
is thus obvious that at least one phase transition occurs 
between these two limiting cases. 

In the presence of the field, A s 's and -Bp's are no longer 
conserved so that H is no longer integrable. However, for 
some special directions of the field, some mappings onto 
well-known problems exist. In the following and without 
loss of generality, we set J = 1/2. 

• h y = - The first simple example is ob- 
tained when the field points in the x (or z) direc- 
tion. In this case, the problem is equivalent to the two- 
dimensional (2D) transverse-field quantum Ising model 
[8, 9] which is known to display a second-order transi- 
tion for h x — 0.1642(2) [ ]. When both x and z compo- 
nents of the field are nonvanishing, the Hamiltonian H is 
equivalent to the 3D classical Z 2 gauge Higgs model [10]. 
In the plane h y = 0, the phase diagram consists of two 
second-order lines which originate from the Ising points 
(h x = and h z = 0) and intersect at a multicritical point 
located at the symmetric point h x — h z = 0.1703(2) [11]. 

• h x = h z = - When the field points in the y di- 
rection, H is self-dual (its spectrum is invariant under 
the exchange h y o J). In addition, it is isospectral to 
the 2D quantum compass model [ ] which is also equiv- 
alent to that of the Xu-Moore model [22]. In this 
first-order transition occurs at the point h y = J [12, 19]. 

Methods : PCUT and iPEPS — Away from these spe- 
cial directions, no mapping onto existing models is known 
so far. To analyze the full phase diagram, we have first 
computed the low-energy spectrum using the PCUT (to- 
gether with the finite-lattice method [23]) in the small- 
field limit, which has already been proven to be very ef- 
ficient in this context [11, 12]. This approach provides a 
natural description in terms of dressed anyonic quasipar- 
ticles in the thermodynamical limit. We focused on the 
ground-state energy per spin eo and the one-quasiparticle 
gap A which have been computed at order 10 and 8, re- 
spectively. The lengthy expressions of these quantities 
can be found in the supplementary material. We em- 
phasize that, at such high orders, eo (A) is determined 



with a relative precision lower than 10 -3 (10 -2 ) for all 
directions of the magnetic field and inside the topological 
phase. Of course, as for any series expansion, such error 
bars can only be roughly estimated using various resum- 
mation schemes (see Ref. 24 for a detailed discussion). 

The PCUT method allows us to determine the set of 
points (h x ,h y ,h z ) where A vanishes and hence where 
there might be a continuous transition. However, we 
know that for h x = h z = 0, the transition is first order 
and thus not detectable by the condition A = 0. This 
is the main reason for using a complementary tool based 
on a variational approach, the so-called iPEPS algorithm, 
which also allows to estimate eo in the thermodynamic 
limit with a rather good accuracy [17-19]. The main pa- 
rameter in this method is the so-called bond dimension 
D of the PEPS tensors [16-18] which drives the amount 
of entanglement of the ansatz states. 

Our main motivation for choosing such ansatz states 
is that eigenstates of the toric code (zero-field limit) are 
described by D — 2 PEPS [25] whereas for J = 0, eigen- 
states of H are D — 1 (completely separable) states. Ob- 
viously, in the large D limit, this variational method gives 
the exact ground state but, in practice, we have checked 
that the difference between D = 2 and D = 3 lies within 
the error bars of the PCUT calculation so that, for the 
sake of simplicity, we restrict our analysis to D = 2 only. 
Once the bond parameter is fixed, one still has the free- 
dom to choose different ansatz states. Here, we choose a 
PEPS structure similar to that proposed in Ref. 17, but 
we allow four different tensors for the four spins of each 
elementary plaquette (instead of two in Ref. ) . Such a 
choice leads to 8D 4 — 1 variational parameters (instead of 
4£> 4 - 1) and thus improves the results. Other technical 
details of the algorithm have also been adapted to tackle 
four-spin interactions. 

One may argue that in order to capture the topological 
properties of the ground state in the general case (such as 
a nontrivial topological entropy [26]), one would need to 
implement some gauge symmetries in the tensor network 
ansatz [27]. But, such properties reflect nonlocal features 
and are not crucial for computing local quantities such 
as the ground-state energy. 

Keeping all these approximations in mind, let us de- 
scribe the general strategy to determine a transition point 
and its nature (first or second order). For a fixed direc- 
tion of the field we wish to compute the critical value 
of the field's strength h beyond which the system is no 
more in a topological phase. To do so, one proceeds in 
three steps : (i) compute the iPEPS ground-state en- 
ergy e PEPS for different values of h by minimizing the 
tensor parameters; (ii) determine the point h* at which 
e iPEPS < e PCUT where e PCUT denotes the PCUT ground- 
state energy; (in) compute the value h c for which the 
one-quasiparticle gap vanishes using the PCUT expres- 
sion of A and resummation techniques. Then two situ- 
ations must be distinguished. Either h* > h c , in which 
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FIG. 1. (Color online) Comparison of iPEPS and PCUT 
ground-state energy for two different field directions. The 
width of the (gray) band denning h c results from differ- 
ent Dlog Pade approximants. Left : h = 0, 1) and 
h* > h c indicating a second-order transition at h c . Right : 
h = h(cos ^,sin |^,cos jg) and h* < h c indicating a first- 
order transition at h* . 



case we can trust the PCUT result and its prediction of 
a second-order transition at h c . The iPEPS approach is 
indeed variational and invalidates the PCUT's prediction 
when eg PEPS < e PCUT . Or h* < h c , in which case a tran- 
sition occurs before the gap A vanishes. This means that 
there are some level crossings due to higher-energy levels 
which are not captured by the PCUT approach, indicat- 
ing a first-order transition confirmed by the discontinuity 
of the slope of the iPEPS energy [see e.g. Fig. 1 (right)]. 
Note that one may indeed directly compute the deriva- 
tive of e PE as a function of h and look for singularities 
but this approach is less precise. Obviously the precision 
in the determination of h* and h c plays a fundamental 
role in this scheme. For a given direction, the maximum 
orders at which we computed eo and A as well as the 
form of the chosen variational states allow us to estimate 
the transition point with an accuracy of a few percent as 
can be seen in Fig. 1. 

Phase diagram — A sketch of the 3D phase diagram 
is shown in Fig. 2 and can be summarized as follows. 
First, we find that the transition point h — (0,1/2,0) 
is part of a 2D first-order transition sheet S±. Second, 
the second-order transition lines of the h y = plane give 
rise to a 2D second-order transition sheet S2 (defined by 
A = 0) when the y-component of the field is nonvanish- 
ing. These sheets that intersect on a nontrivial line define 
the boundaries of the topological phase. Given the diffi- 
culty for investigating the full 3D space with iPEPS, we 
focused on some special planes in which we determined 
the coordinates of the intersection point of S\ and S2. 
For instance, in the (0,h y ,h z ) plane, we found that this 
intersection occurs around the point h — (0,0.49,0.11). 
When the transition is second order, the gap is expected 
to behave as A ~ (h — h c ) ZL ' in the vicinity of the critical 



FIG. 2. (Color online) Sketch of the 3D phase diagram. Dots 
correspond to Ising points and the diamond is the self-dual 
point of the h y line. Green lines are the intersections of the 
first-order sheet Si and the second-order sheet S2 (computed 
from the bare series given in supplementary material). The 
multicritical line h x — h z with continuously varying critical 
exponents is shown as a thick (red) line. 

point h c . Note that here, we do not have access to the 
dynamical exponent z and to the correlation length expo- 
nent v independently but only to their product. For all 
investigated directions, we found that zv was compatible 
with the well-established Ising value zv — 0.630(1). This 
leads us to conclude that 1S2 lies in the Ising universality 
class (as was already found in the plane h y = [10, 11]) 
for all directions except for the special case h x = h z . 

The multicritical line — As discussed in [10, 11] for 
h y = 0, the two second-order transition lines merge in a 
multicritical point at h x = h z for which the gap exponent 
is clearly different from the Ising value. The most impor- 
tant result of the present study is that when h y 7^ 0, 
this multicritical point gives rise to a multicritical line 
on which this exponent varies continuously. First of all, 
let us point out that the multicritical line intersects 1S1 
around the point h — (0.17, 0.46, 0.17). Once again these 
values are obtained with a relative precision of a few per- 
cent. Along this multicritical line, we have computed 
the exponent zv using standard resummation techniques 
based on Dlog Pade approximants (see Ref. 24 for de- 
tails). Our results are displayed in Fig. 3 and show that 
this exponent varies from 0.69 at h y = up to a value 
close to 1 at h y — 0.46 along this line. Except in the range 
h y G [0.20, 0.35], one gets a rather good convergence sug- 
gesting that divergencies observed in this region are due 
to spurious poles in the Dlog Pade approximants. We 
thus conjecture that zv varies continuously and that its 
variation of ~ 50% cannot be attributed to extrapola- 
tion errors and reveals a new universality class. Since it 
is not associated to a symmetry breaking but rather re- 
flects the breakdown of a topological phase, we will call 
it topological. 

At this stage, it is difficult to determine the key in- 
gredients for a system to belong to this class (since we 
do not have any local order parameter) but it is likely 
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FIG. 3. (Color online) Critical exponent zv as a function of 
h y along the line h x = h z computed for various Dlog Pade ap- 
proximants [m, n]. Strange behaviors near h y ~ 0.3 are likely 
due to spurious pole structures and should not be considered 
as relevant. 



that the mutual semionic statistics of charges and fluxes 
is one of them. More generally, let us underline that con- 
tinuously varying critical exponents are not common in 
two-dimensional quantum systems. During the comple- 
tion of this work, some conformal quantum critical lines 
in 2+1 dimensions have been proposed [28, 29] but their 
relevance for the toric code in a magnetic field is still an 
open question. 

Discussion and outlook — In the present work, we 
have determined the boundaries of the topological phase 
of the toric code in a field using two state-of-the-art and 
complementary methods. This topological "bubble" is 
made of first-order and second-order sheets. Interest- 
ingly, second-order transitions seem to be in the Ising 
universality class except on a multicritical line on which 
the gap vanishes with continuously varying exponents 
giving rise to a new "topological" universality class. Of 
course, it would also be valuable to study the large-field 
limit of this model to investigate the outer part of 
the bubble. Notably the fate of the first-order line 
observed in the h y = plane [10, 11] is an interesting 
question. Finally, a complete understanding of the 
low-energy spectrum of the topological phase certainly 
requires the study of bound states as already seen in the 
transverse-field case [12]. 
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SUPPLEMENTARY MATERIAL 

Here are the series expansions obtained using the PCUT method in the small-field limit h x ,h y ,h z <C J. Setting 
Sk = h x + h z , P 2 k — h x h z and J — 1/2, the ground-state energy per spin eo at order 10 reads 

1 S 2 h* X5Sj 7S 2 hl P 4 _ I3h 4 y _ U7S 6 _ 3716^ 113S , 2 P 4 _ lOA5S 2 h y 2003P 4 h 2 y 
C ° ~ ~2 ~ Y ~ T 8 32~ + T ~ 192 ~ ~~ 8 128 + 32 3456 + 384 

Wh y 1800358 19548795 6 /i2 66855 4 P 4 340541755 4 ^ 1468615 , 2 P 4 ^ 153435495 2 ^ 
~ 3072 ~ 64 ~ 36864 + 128 3981312 + ~ 2304 26542080 

20869P 8 5020085P 4 /^ _ 163885^ _ 5420775510 _ 15634595235 8 ftg 39524033S , 6 P 4 
+ 384 + 497664 1769472 1024 1327104 + 36864 

111510540942756^ 100582354455 4 P 4 /i2 42196408354975 4 ^ 56509255 2 P 8 2085409756352 P±h\ 
5733089280 + " 7962624 ~ 191102976000 _ 6912 " 143327232 



(1) 



4838909402815 2 /i® 1202498305P 8 /j2 1994817656221P 4 ^ 186734746441/iJ 
382205952000 + 1990656 + 71663616000 1146617856000 ' 

Similarly, for < h x < h z , and the one-quasiparticle (dressed charge) gap A at order 8 reads 

11 1 K 07 

A = 1 — ih z - hy - 4h 2 z + 2h 2 x h z + —h 2 y h z - 12h 3 z + hh x + 17h 2 x h 2 y - —h y + 3h 2 h 2 z - 9h 2 y h 2 z - 36/i* + —h%h z 

+ X lh% + \h 2 x h 2 y h z + ^h 2 x hl + ^h 2 y hl - 178*« + 92/4 + ^Khi + nh 2 z hi + ^§Kh 2 x + GShft 
1305 2 2 2 575 . 2625 . 7971 , a , 4 135619 4 2 495 . 1142149 4 2 3031 a 4 

- 384^ " ^ " -aT^i " + T^*' + "W^*' - 13824^ 

799973 925 13807 1782929,4,3 + 28633 _ 238621 _ 14771 35649 

110592 y z 4 x z 48 x v z 20736 v z 64 z 2 1152 y 2 4 2 16 2 

7715431 , B , 2 3032191, 4 , 4 98263727 , 2 , fi 26492351 , R 80999, fil2 2199571 , 4 , 2l 2 

H /i^/i?, H H hlht hi H hlh 2 + fti/i?,^ 

3072 x v 31104 * y 3981312 x v 7962624 » 96 x 2 4608 2 y 2 

24547709,2, 4 ,2 1495320677 , R , , 19263, 4 , 4 5186533 , 2 , 2 , 4 1760584999 , 4 , 4 118029, 2lfi 

H hihth, Khi H H hihih; h M/i, 

165888 x y 2 19906560 y 2 16 x 2 1728 2 a 2 1990656 2/2 64 

4663837 , . « 940739 , 



1728 a 2 64 1 7 

The (dressed flux) gap for h x > h z is straightforwardly obtained by exchanging h x and h z in this expression. 



Errata : 

- For h y = 0, one recovers expressions given in Eq. (8) of Ref. 11 up to a typo : the term proportional to (h x + h s z ) 
is missing. 

- For h x = h z = 0, one recovers expressions given in Eq. (4) of Ref. 12 up to a typo : the term proportional to t w 
must be corrected by a factor 2. 



